/*  partition.cpp
 * Copyright (C) 2010, Francisco Claude & Patrick Nicholson, all rights reserved.
 *
 * Francisco Claude & Patrick Nicholson <(fclaude|p3nichol)@cs.uwaterloo.ca>
 *
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include <iostream>
#include <iomanip>
//#include <stdio.h>
#include <string>
#include <set>

#include <vector>
#include <map>
#include <list>
#include <limits>
#include <cmath>
#include <cstdlib>
#include <fstream>

using namespace std;

/**
 * A simple class for storing 2D points.  Complete with overloaded
 * comparator and stream output.
 */
class Point2 {
public:

  double x, y;
  Point2(const Point2 &p) { x = p.x; y = p.y; }
  Point2(double _x, double _y) { x = _x; y = _y; }
  Point2() { x = 0; y = 0; }

  bool operator <(const Point2& b) const {
    return x < b.x || (x == b.x && y < b.y);
  }

  friend ostream& operator<< (ostream &out, Point2 &cPoint);
};

ostream& operator<< (ostream &out, Point2 &p) {
  out << "(" << p.x << "," << p.y << ")";
  return out;
}


/**
 * Layer class, used to encapsulate the creation of layer structures
 * on the input data set, and also allow the extraction of the longest
 * increasing subsequence from the data.  Currently only decending
 * layers are supported.
 */
class Layers {

private:

  int n;
  vector<list<Point2> > layerlist;	

public:

  static const int ASCENDING = 0;
  static const int DESCENDING = 1;
  int direction;
 
  /**
   * The Layer constructor takes a direction and an STL set of points
   * as input, and then builds the set of layers in the specified
   * direction using the algorithm described in Yehuda and Fogel's
   * paper.  Assumes thePoints is non-empty.
   */
  Layers(int dir, set<Point2> &thePoints) {
    
    direction = dir;

    Point2 *P = new Point2[thePoints.size()];
    n = 0;

    for (set<Point2>::iterator sit = thePoints.begin(); sit != thePoints.end(); sit++) {
      P[n++] = *sit;
    }

    double *A = new double[n+1];
    A[0] = numeric_limits<double>::max();

    int l = 0;
    for (int i = n-1; i >= 0; i--) {
      int j = 0;
      int j_prime = 0;
      for(; j <= l; j++) {				
	if ((direction == DESCENDING && A[j] > P[i].y ) || 
	    (direction == ASCENDING && A[j] < P[i].y)) {
	  j_prime = j;
	}
      }
      j = j_prime;

      if (j == l) {
	layerlist.push_back(list<Point2>());
	l++;
      }

      list<Point2> *L = &layerlist[j];
      L->push_front(Point2(P[i]));
      A[j+1] = P[i].y;
    }
    delete [] P;
    delete [] A;
  }

  /**
   * Delete the points in the layer structure.
   */
  ~Layers() {
    layerlist.empty();
  }

  /**
   * Return the number of layers.
   */
  int size() {
    return layerlist.size();
  }

  /**
   * Identify and return a list of points that describe an antichain
   * (that is a set of points going in the opposite direction of the
   * layers) that has length k.  Obviously k cannot be greater than
   * the number of layers or this will not work.
   */
  list<Point2> getIncreasingSubsequence(int k) {

    list<Point2> lis;

    if (layerlist.size() > 0) {
      int l = layerlist.size() - 1;
      Point2 q_prime = layerlist[l].front();
      lis.push_front(Point2(q_prime));
      for (int i = l - 1; i >= (l - k + 1) && l >= 0; i--) {
	for (list<Point2>::iterator lit = layerlist[i].begin(); lit != layerlist[i].end(); lit++) {
	  Point2 q = *lit;
	  if (q.x >= q_prime.x && 
	      ((direction == DESCENDING && q.y >= q_prime.y) ||
	       (direction == ASCENDING && q.y <= q_prime.y))) {
	    lis.push_front(Point2(q));
	    q_prime = Point2(q);
	    break;
	  }
	}
      }
    }
    return lis;
  }
};


/**
 * Partition the set "points" into two sets, "dscPoints" and
 * "ascPoints", using (a not-so-efficient-version of) the algorithm
 * described in Yang et al. 2007 "A Comparative Sutdy of Efficient
 * Algorithms for Partitioning a Sequence into Monotone Subsequences."
 * It is less efficient since instead of using the layer update method
 * described by Yehuda and Fogel, it just recomputed the layers at
 * each step in the algorithm.  Thus, the complexity is
 * O(n*lg(n)*sqrt(n)), instead of O(n*sqrt(n)).
 */
void partitionBing(set<Point2> &points, set<Point2> &dscPoints, set<Point2> &ascPoints) {
		
  if (points.empty()) {
    return;
  }
   
  dscPoints.clear();
  ascPoints.clear();

  for (set<Point2>::iterator sit = points.begin(); sit != points.end(); sit++) {
    dscPoints.insert(*sit);
  }

  
  int startCount = -1;
  int goal = -1;
  
  int sequences = 0;
  int K = (int) (sqrt(2.0 * points.size() + 0.25) - 0.5);
   
  cerr << setw(3) << 0 << "%";

  for (int k = K; k > 0; k--) {

    Layers dscLayers = Layers(Layers::DESCENDING, dscPoints);
    int l = dscLayers.size();

    if (startCount == -1) {
      startCount = l;
      goal = startCount - K;
    }
    int progress = (int) ((((float)startCount - (l + sequences)) / goal) * 100.0);
    cerr << "\b\b\b\b" << setw(3) << progress << "%";
   
    if (l <= k) {
      sequences += l;
      break;				
    } else {
      list<Point2> lis = dscLayers.getIncreasingSubsequence(k);
      for (list<Point2>::iterator it = lis.begin(); it != lis.end(); it++) {
	Point2 p = *it;
	ascPoints.insert(p);
	dscPoints.erase(p);
      }
      sequences++;
      if (dscPoints.size() == 0) {
	break;
      }				
    }
  }

  cerr <<  endl;

  cerr << "Number of points:    " << points.size() << endl;
  cerr << "Number of sequences: " << sequences << endl;
  cerr << "Bound on sequences:  " << K << endl;
}


/**
 * Partition the set of points into two sets: dscPoints which are
 * points that are in descending chains, and ascPoints which are
 * points that are in ascending chains.  The algorithm is described by
 * Fomin, Kratsch, and Novelli in their paper "Approximating minimum
 * cocolorings.  In the paper version they greedly select the largest
 * clique and independent sets from a graph to obtain a cocoloring
 * with size at most log(n) * OPT.  Here, we can simply greedily
 * select maximum chains and anti-chains until no points are left.
 */
void partitionGreedy(set<Point2> &points, set<Point2> &dscPoints, set<Point2> &ascPoints) {
		
  if (points.empty()) {
    return;
  }
  
  set<Point2> pointsCopy;
  pointsCopy.insert(points.begin(), points.end());

  dscPoints.clear();
  ascPoints.clear();

  int sequences = 0;

  int startPoints = pointsCopy.size();

  cerr << setw(3) << 0 << "%";

  while (true) {
    cerr << "\b\b\b\b" << setw(3) << 
      (int) (100 * (1.0 - ((1.0 * pointsCopy.size()) / startPoints))) << "%";

    Layers dscLayers = Layers(Layers::DESCENDING, pointsCopy);
    Layers ascLayers = Layers(Layers::ASCENDING, pointsCopy);

    list<Point2> lis;
    if (ascLayers.size() > dscLayers.size()) {
      lis = ascLayers.getIncreasingSubsequence(ascLayers.size());
      for (list<Point2>::iterator it = lis.begin(); it != lis.end(); it++) {
	dscPoints.insert(*it);
	pointsCopy.erase(*it);
      }
    } else {
      lis = dscLayers.getIncreasingSubsequence(dscLayers.size());
      for (list<Point2>::iterator it = lis.begin(); it != lis.end(); it++) {
	ascPoints.insert(*it);
	pointsCopy.erase(*it);
      }
    }
    sequences++;

    if (pointsCopy.empty()) {
      break;
    }
  }
  cerr << "\b\b\b\b" << "100%" << endl;
  cerr << "Number of points:    " << points.size() << endl;
  cerr << "Number of sequences: " << sequences << endl;

}


/**
 * Read the input file, which contains lines of the format:
 *
 * <point number> <x-coordinate> <y-coordinate>
 *
 * where the x and y coordinates are floating point numbers, and dump
 * the points into an STL set. Returns true if the file as read
 * without errors.
 *
 */
bool loadPoints(string inFilePath, set<Point2> &points) {
  ifstream infile (inFilePath.c_str());
  string line;
  int lineno = 1;
  if (infile.is_open()) {
    while (! infile.eof() ) {
      getline (infile,line);
      double px, py;
      int match = sscanf(line.c_str(), "%*d %lf %lf", &px, &py);
      if (match == 2) {
	points.insert(Point2(px,py));
      } else {
	//cerr << "Error parsing file on line " << lineno << ": " << line << endl;
	//cerr << "Only matched " << match << " tokens. Skipping line." << endl;
      }
      lineno++;
    }
    infile.close();
    return true;
  }
  return false;
}


/**
 * Write the points from dscPoints and ascPoints into a file using
 * the same format as the input file, except adding a second column:
 *
 * <point number> <ascending/descending> <x-coordinate> <y-coordinate>
 *
 * Returns true if the file is created and written to without errors.
 */
bool writePartitionedPoints(string outFilePath, set<Point2> &dscPoints, set<Point2> &ascPoints) {

  ofstream outfile(outFilePath.c_str());
  outfile.precision(14);
  outfile << dscPoints.size() << " "  << ascPoints.size() << endl;
  int pNum = 0;
  if (outfile.is_open()) {
    for (set<Point2>::iterator it = dscPoints.begin(); it != dscPoints.end(); it++) {
      Point2 p = *it;
      outfile << pNum++ << " D " << p.x << " " << p.y << endl;
    }

    for (set<Point2>::iterator it = ascPoints.begin(); it != ascPoints.end(); it++) {
      Point2 p = *it;
      outfile << pNum++ << " A " << p.x << " " << p.y << endl;
    }
    outfile.close();
    return true;
  } else {
    return false;
  }
}


#define GREEDY 0
#define BING_ET_AL 1

/**
 * Main function: read the input file of points, partition them using
 * the computePartition function, then write the output to the
 * specified output file.
 */
int main(int argc, char **argv) {

  set<Point2> points, dscPoints, ascPoints;
  
  int method = GREEDY;
  char* progName = argv[0];
  argc--;
  argv++;

  if(argc==0) {
    cerr << "Usage: " << progName << " <infile> <outfile>" << endl;
    exit(1);
  }
  
  /*if (string(argv[0]) == "-greedy") {
    method = GREEDY;
    argc--;
    argv++;
  } else if (string(argv[0]) == "-bing") {
    method = BING_ET_AL;
    argc--;
    argv++;
  }*/

  if (argc != 2) {
    cerr << "Usage: " << progName << " <infile> <outfile>" << endl;
    exit(1);
  }

  string inFilePath(argv[0]);
  string outFilePath(argv[1]);

  if (loadPoints(inFilePath, points)) {
    cerr << "Opened file \"" << inFilePath << "\"" << endl;
    switch(method) {
    case BING_ET_AL: 
      partitionBing(points, dscPoints, ascPoints);
      break;
    case GREEDY:
      partitionGreedy(points, dscPoints, ascPoints);
      break;
    default:
      cerr << "Selected partitioning method does not exist." << endl;
      exit(1);
    }
    if (writePartitionedPoints(outFilePath, dscPoints, ascPoints)) {
      cerr << "Wrote output to file \"" << outFilePath << "\"" << endl;
    } else {
      cerr << "Failed to write output to file \"" << outFilePath  << "\"" << endl;
      return 1;
    }
  } else {
    cerr << "Error opening file \"" << inFilePath << "\" for reading." << endl;
    return 1;
  }
  return 0;
}

